Pylops - Marchenko redatuming

Author: M.Ravasi

In this notebook I will show you how to use the MDC and Marchenko linear operators of Pylops to perform Marchenko redatuming by inversion - old implementation, used to check backcompatibility with versions <v1.5.0, head over to Marchenko for a more complete set of examples.

In [1]:
%load_ext memory_profiler
%load_ext autoreload
%autoreload 2
%matplotlib inline

import warnings
warnings.filterwarnings('ignore')

import numpy as np
import matplotlib.pyplot as plt
import scipy as sp

from scipy.sparse import csr_matrix, vstack
from scipy.linalg import lstsq, solve
from scipy.sparse.linalg import LinearOperator, cg, lsqr
from scipy.signal import convolve, filtfilt

from pylops.utils                      import dottest
from pylops.utils.wavelets             import *
from pylops.utils.seismicevents        import *
from pylops.utils.tapers               import *
from pylops.basicoperators             import *
from pylops.signalprocessing           import *
from pylops.waveeqprocessing.mdd       import *
from pylops.waveeqprocessing.marchenko import *
from pylops.optimization.leastsquares  import *

# just for comparison
from marchenko_old import Marchenko as Marchenko

Inputs

Input parameters

In [2]:
inputfile = '../../../pylops/testdata/marchenko/input.npz' # choose file in testdata folder of repo

vel = 2400.0        # velocity
toff = 0.045        # direct arrival time shift
nsmooth = 10        # time window smoothing 
nfmax = 500         # max frequency for MDC (#samples)
nstaper = 11        # source/receiver taper lenght
n_iter = 15         # iterations

jr = 1              # subsampling in r
js = 1              # subsampling in s

Load input

In [3]:
inputdata = np.load(inputfile)

Read and visualize geometry

In [4]:
# Receivers
r = inputdata['r'][:,::jr]
nr = r.shape[1]
dr = r[0,1]-r[0,0]

# Sources
s = inputdata['s'][:,::js]
ns = s.shape[1]
ds = s[0,1]-s[0,0]

# Virtual points
vs = inputdata['vs']

# Density model
rho = inputdata['rho']
z, x = inputdata['z'], inputdata['x']

plt.figure(figsize=(18,9))
plt.imshow(rho, cmap='gray', extent = (x[0], x[-1], z[-1], z[0]))
plt.scatter(s[0, 5::10], s[1, 5::10], marker='*', s=150, c='r', edgecolors='k')
plt.scatter(r[0, ::10],  r[1, ::10], marker='v', s=150, c='b', edgecolors='k')
plt.scatter(vs[0], vs[1], marker='.', s=250, c='m', edgecolors='k')
plt.axis('tight')
plt.xlabel('x [m]'),plt.ylabel('y [m]'),plt.title('Model and Geometry')
plt.xlim(x[0], x[-1]);

Read data

In [5]:
# time axis
t = inputdata['t']
ot, dt, nt = t[0], t[1]-t[0], len(t)

# data
R = inputdata['R'][::js, ::jr]
R = np.swapaxes(R, 0, 1)

# tapering
taper = taper3d(nt, [ns, nr], [nstaper, nstaper], tapertype='hanning')
R = R*taper

Read subsurface fields and wavelet to apply to subsurface fields

In [6]:
Gsub = inputdata['Gsub'][:, ::jr]
G0sub = inputdata['G0sub'][:, ::jr]
wav = inputdata['wav']
wav_c = np.argmax(wav)

# convolve with wavelet
Gsub = np.apply_along_axis(convolve, 0, Gsub, wav, mode='full')
Gsub = Gsub[wav_c:][:nt]
G0sub = np.apply_along_axis(convolve, 0, G0sub, wav, mode='full') 
G0sub = G0sub[wav_c:][:nt]
In [7]:
fig, axs = plt.subplots(1, 3, sharey=True, figsize=(15, 9))
axs[0].imshow(R[20].T, cmap='gray', vmin=-1e-2, vmax=1e-2, extent=(r[0,0], r[0,-1], t[-1], t[0]))
axs[0].set_title('R shot=0'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(1.5, 0)
axs[1].imshow(R[ns//2].T, cmap='gray', vmin=-1e-2, vmax=1e-2, extent=(r[0,0], r[0,-1], t[-1], t[0]))
axs[1].set_title('R shot=%d' %(ns//2)), axs[1].set_xlabel(r'$x_R$')
axs[1].axis('tight')
axs[1].set_ylim(1.5, 0)
axs[2].imshow(R[-20].T, cmap='gray', vmin=-1e-2, vmax=1e-2, extent=(r[0,0], r[0,-1], t[-1], t[0]))
axs[2].set_title('R shot=%d' %ns), axs[2].set_xlabel(r'$x_R$')
axs[2].axis('tight');
axs[2].set_ylim(1.5, 0)

fig, axs = plt.subplots(1, 2, sharey=True, figsize=(12, 9))
axs[0].imshow(Gsub, cmap='gray', vmin=-1e5, vmax=1e5, extent=(r[0,0], r[0,-1], t[-1], t[0]))
axs[0].set_title('G'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(1.5, 0)
axs[1].imshow(G0sub, cmap='gray', vmin=-1e5, vmax=1e5, extent=(r[0,0], r[0,-1], t[-1], t[0]))
axs[1].set_title('G0'), axs[1].set_xlabel(r'$x_R$'), axs[1].set_ylabel(r'$t$')
axs[1].axis('tight')
axs[1].set_ylim(1.5, 0);

Marchenko preparation

Create window

In [8]:
# direct arrival window - traveltime
directVS = np.sqrt((vs[0]-r[0])**2+(vs[1]-r[1])**2)/vel
directVS_off = directVS - toff

# window
idirectVS_off = np.round(directVS_off/dt).astype(np.int)
w = np.zeros((nr, nt))
for ir in range(nr):
    w[ir, :idirectVS_off[ir]]=1            
w = np.hstack((np.fliplr(w), w[:, 1:]))

if nsmooth>0:
    smooth=np.ones(nsmooth)/nsmooth
    w  = filtfilt(smooth, 1, w)    
    
fig, ax = plt.subplots(1, 1,  sharey=True, figsize=(5, 9))
im = ax.imshow(w.T, cmap='gray', extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
ax.plot(r[0], directVS,'b'),ax.plot(r[0], -directVS,'b')
ax.set_title('Window'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
ax.axis('tight')
fig.colorbar(im, ax=ax);

Create analytical direct wave

In [9]:
G0sub_ana = directwave(wav, directVS, nt, dt, nfft=2**11)

fig = plt.figure(figsize=(15,9))
ax1 = plt.subplot2grid((1, 5), (0, 0), colspan=2)
ax2 = plt.subplot2grid((1, 5), (0, 2), colspan=2)
ax3 = plt.subplot2grid((1, 5), (0, 4))

ax1.imshow(G0sub/G0sub.max(), cmap='gray', vmin=-1, vmax=1, 
           extent=(r[0,0], r[0,-1], t[-1], t[0]))
ax1.set_title(r'$G0_{FD}$'), ax1.set_xlabel(r'$x_R$'), ax1.set_ylabel(r'$t$')
ax1.axis('tight')
ax2.imshow(G0sub_ana/G0sub_ana.max(), cmap='gray', vmin=-1, vmax=1, 
           extent=(r[0,0], r[0,-1], t[-1], t[0]))
ax2.set_title(r'$G0_{Ana}$'), ax2.set_xlabel(r'$x_R$'), ax2.set_ylabel(r'$t$')
ax2.axis('tight')

ax3.plot(G0sub[:, nr//2]/G0sub.max(), t, 'r', lw=5)
ax3.plot(G0sub_ana[:, nr//2]/G0sub_ana.max(), t, 'k', lw=3)
ax3.set_ylim(1.2, 0);

For now we will use the direct wave obtained via finite-difference, but we will show later that Marchenko redatuming can converge to the correct solution also when using an analytical direct wave.

Inversion

In [10]:
# Add negative time to operators
Rtwosided = np.concatenate((np.zeros((nr, ns, nt-1)), R), axis=-1)
R1twosided = np.concatenate((np.flip(R, axis=-1), 
                            np.zeros((nr, ns, nt-1))), axis=-1)

Rtwosided_fft = np.fft.rfft(Rtwosided, 2*nt-1, axis=-1)/np.sqrt(2*nt-1)
Rtwosided_fft = Rtwosided_fft[...,:nfmax]
R1twosided_fft = np.fft.rfft(R1twosided, 2*nt-1, axis=-1)/np.sqrt(2*nt-1)
R1twosided_fft = R1twosided_fft[...,:nfmax]

# Operators
Rop = MDC(Rtwosided_fft, nt=2*nt-1, nv=1, dt=dt, dr=ds, twosided=True, dtype='complex64')
dottest(Rop, (2*nt-1)*nr, (2*nt-1)*ns, verb=True)

R1op = MDC(R1twosided_fft, nt=2*nt-1, nv=1, dt=dt, dr=ds, twosided=True, dtype='complex64')
dottest(R1op, (2*nt-1)*nr, (2*nt-1)*ns, verb=True)

# Input focusing function
fd_plus =  np.concatenate((np.fliplr(G0sub.T), np.zeros((nr, nt-1))), axis=-1)

Dot test passed, v^T(Opu)=-86.819385 - u^T(Op^Tv)=-86.819385
Dot test passed, v^T(Opu)=57.118539 - u^T(Op^Tv)=57.118539

Create Marchenko operator

In [11]:
Wop = Diagonal(w.flatten())
Iop = Identity(nr*(2*nt-1))

Mop = VStack([HStack([Iop, -1*Wop*Rop]),
              HStack([-1*Wop*R1op, Iop])])*BlockDiag([Wop, Wop])
Gop = VStack([HStack([Iop, -1*Rop]),
              HStack([-1*R1op, Iop])])

dottest(Gop, 2*nr*(2*nt-1), 2*nr*(2*nt-1), verb=True)
dottest(Mop, 2*nr*(2*nt-1), 2*nr*(2*nt-1), verb=True);

Dot test passed, v^T(Opu)=-112.217832 - u^T(Op^Tv)=-112.217832
Dot test passed, v^T(Opu)=-87.473918 - u^T(Op^Tv)=-87.473918

Run standard redatuming as benchmark

In [12]:
p0_minus = Rop * fd_plus.flatten()
p0_minus = p0_minus.reshape(nr, (2*nt-1))

Create data, adjoint and inverse focusing functions

In [13]:
d = Wop*Rop*fd_plus.flatten()
d = np.concatenate((d.reshape(nr, 2*nt-1), np.zeros((nr, 2*nt-1))))

f1_adj = Mop.H*d.flatten()
f1_inv = lsqr(Mop, d.flatten(), iter_lim=n_iter, show=True)[0]

f1_adj = f1_adj.reshape(2*nr, (2*nt-1))
f1_inv = f1_inv.reshape(2*nr, (2*nt-1))

 
LSQR            Least-squares solution of  Ax = b
The matrix A has   322998 rows  and   322998 cols
damp = 0.00000000000000e+00   calc_var =        0
atol = 1.00e-08                 conlim = 1.00e+08
btol = 1.00e-08               iter_lim =       15
 
   Itn      x[0]       r1norm     r2norm   Compatible    LS      Norm A   Cond A
     0  0.00000e+00   2.983e+07  2.983e+07    1.0e+00  3.5e-08
     1  0.00000e+00   1.311e+07  1.311e+07    4.4e-01  9.2e-01   1.1e+00  1.0e+00
     2  0.00000e+00   7.406e+06  7.406e+06    2.5e-01  3.9e-01   1.8e+00  2.2e+00
     3  0.00000e+00   5.479e+06  5.479e+06    1.8e-01  3.3e-01   2.1e+00  3.4e+00
     4  0.00000e+00   3.659e+06  3.659e+06    1.2e-01  3.3e-01   2.5e+00  5.2e+00
     5  0.00000e+00   2.780e+06  2.780e+06    9.3e-02  2.6e-01   2.9e+00  6.8e+00
     6  0.00000e+00   2.244e+06  2.244e+06    7.5e-02  2.3e-01   3.2e+00  8.5e+00
     7  0.00000e+00   1.498e+06  1.498e+06    5.0e-02  2.5e-01   3.6e+00  1.1e+01
     8  0.00000e+00   1.105e+06  1.105e+06    3.7e-02  1.9e-01   3.9e+00  1.3e+01
     9  0.00000e+00   8.988e+05  8.988e+05    3.0e-02  1.6e-01   4.2e+00  1.4e+01
    10  0.00000e+00   6.714e+05  6.714e+05    2.3e-02  1.7e-01   4.4e+00  1.6e+01
    11  0.00000e+00   5.040e+05  5.040e+05    1.7e-02  1.5e-01   4.6e+00  1.8e+01
    12  0.00000e+00   4.189e+05  4.189e+05    1.4e-02  1.4e-01   4.8e+00  2.0e+01
    13  0.00000e+00   3.357e+05  3.357e+05    1.1e-02  1.4e-01   5.0e+00  2.2e+01
    14  0.00000e+00   2.718e+05  2.718e+05    9.1e-03  1.1e-01   5.3e+00  2.4e+01
    15  0.00000e+00   2.424e+05  2.424e+05    8.1e-03  7.9e-02   5.5e+00  2.6e+01
 
LSQR finished
The iteration limit has been reached                      
 
istop =       7   r1norm = 2.4e+05   anorm = 5.5e+00   arnorm = 1.0e+05
itn   =      15   r2norm = 2.4e+05   acond = 2.6e+01   xnorm  = 3.5e+07

Add initial guess to estimated focusing functions

In [14]:
f1_adj_tot = f1_adj + np.concatenate((np.zeros((nr, 2*nt-1)),
                                      np.concatenate((np.fliplr(G0sub.T), 
                                                      np.zeros((nr, nt-1))), axis=-1)))

f1_inv_tot = f1_inv + np.concatenate((np.zeros((nr, 2*nt-1)),
                                      fd_plus))

Estimate Green's functions

In [15]:
g_adj = Gop*f1_adj_tot.flatten()
g_inv = Gop*f1_inv_tot.flatten()

g_adj = g_adj.reshape(2*nr, (2*nt-1))
g_inv = g_inv.reshape(2*nr, (2*nt-1))

Extract up and down focusing and Green's functions from model vectors

In [16]:
f1_adj_minus, f1_adj_plus =  f1_adj_tot[:nr], f1_adj_tot[nr:]
f1_inv_minus, f1_inv_plus =  f1_inv_tot[:nr], f1_inv_tot[nr:]

g_adj_minus, g_adj_plus =  -g_adj[:nr], np.fliplr(g_adj[nr:])
g_inv_minus, g_inv_plus =  -g_inv[:nr], np.fliplr(g_inv[nr:])

g_inv_tot = g_inv_minus + g_inv_plus

Visualization

In [17]:
fig, axs = plt.subplots(1, 5, sharey=True, figsize=(16, 7))
axs[0].imshow(d.T, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[0].set_title('d'), axs[0].set_xlabel(r'$x_R$'), axs[3].set_ylabel(r'$t$')
axs[0].plot(r[0], directVS,'b'),axs[0].plot(r[0], -directVS,'b')
axs[0].plot(r[0]+r[0,-1], directVS,'b'),axs[0].plot(r[0]+r[0,-1], -directVS,'b')
axs[0].axis('tight')
axs[0].set_ylim(1, -1)
axs[1].imshow(f1_adj_tot.T, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$f_{adj}$'), axs[0].set_xlabel(r'$x_R$')
axs[1].plot(r[0], directVS,'b'),axs[1].plot(r[0], -directVS,'b')
axs[1].plot(r[0]+r[0,-1], directVS,'b'),axs[1].plot(r[0]+r[0,-1], -directVS,'b')
axs[1].axis('tight')
axs[1].set_ylim(1, -1);
axs[2].imshow(f1_inv_tot.T, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$f_{inv}$'), axs[0].set_xlabel(r'$x_R$')
axs[2].plot(r[0], directVS,'b'),axs[2].plot(r[0], -directVS,'b')
axs[2].plot(r[0]+r[0,-1], directVS,'b'),axs[2].plot(r[0]+r[0,-1], -directVS,'b')
axs[2].axis('tight')
axs[2].set_ylim(1, -1);
axs[3].imshow(g_adj.T, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[3].set_title(r'$g_{adj}$'), axs[0].set_xlabel(r'$x_R$')
axs[3].plot(r[0], directVS,'b'),axs[3].plot(r[0], -directVS,'b')
axs[3].plot(r[0]+r[0,-1], directVS,'b'),axs[3].plot(r[0]+r[0,-1], -directVS,'b')
axs[3].axis('tight')
axs[3].set_ylim(1, -1);
axs[4].imshow(g_inv.T, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], 2*r[0,-1], t[-1], -t[-1]))
axs[4].set_title(r'$g_{inv}$'), axs[0].set_xlabel(r'$x_R$')
axs[4].plot(r[0], directVS,'b'),axs[4].plot(r[0], -directVS,'b')
axs[4].plot(r[0]+r[0,-1], directVS,'b'),axs[4].plot(r[0]+r[0,-1], -directVS,'b')
axs[4].axis('tight')
axs[4].set_ylim(1, -1);
In [18]:
fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_minus.T, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].plot(r[0], directVS,'b'),axs[0].plot(r[0], -directVS,'b')
axs[0].axis('tight')
axs[0].set_ylim(1, -1);
axs[1].imshow(f1_inv_minus.T, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$f^-$'), axs[0].set_xlabel(r'$x_R$')
axs[1].plot(r[0], directVS,'b'),axs[1].plot(r[0], -directVS,'b')
axs[1].axis('tight')
axs[1].set_ylim(1, -1);
axs[2].imshow(f1_inv_plus.T, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$f^+$'), axs[0].set_xlabel(r'$x_R$')
axs[2].plot(r[0], directVS,'b'),axs[2].plot(r[0], -directVS,'b')
axs[2].axis('tight')
axs[2].set_ylim(1, -1);
In [19]:
fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_minus.T, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(2, 0);
axs[1].imshow(g_inv_minus.T, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$g^-$'), axs[0].set_xlabel(r'$x_R$')
axs[1].axis('tight')
axs[1].set_ylim(2, 0);
axs[2].imshow(g_inv_plus.T, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$g^+$'), axs[0].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(2, 0);
In [20]:
fig = plt.figure(figsize=(15,9))
ax1 = plt.subplot2grid((1, 5), (0, 0), colspan=2)
ax2 = plt.subplot2grid((1, 5), (0, 2), colspan=2)
ax3 = plt.subplot2grid((1, 5), (0, 4))

ax1.imshow(Gsub, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], t[0]))
ax1.set_title(r'$G_{true}$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
ax1.axis('tight')
ax1.set_ylim(1.2, 0)
ax2.imshow(g_inv_tot.T, cmap='gray', vmin=-5e5, vmax=5e5, 
              extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
ax2.set_title(r'$G_{est}$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
ax2.axis('tight')
ax2.set_ylim(1.2, 0)

ax3.plot(Gsub[:, nr//2]/Gsub.max(), t, 'r', lw=5)
ax3.plot(g_inv_tot[nr//2, nt-1:]/g_inv_tot.max(), t, 'k', lw=3)
ax3.set_ylim(1.2, 0);

Examples of Marchenko class usage for single virtual point

Let's start by feeding the reflection response directly in the freqency domain

In [21]:
MarchenkoWM = Marchenko(Rtwosided_fft, R1=R1twosided_fft, nt=nt, dt=dt, dr=dr, 
                        nfmax=nfmax, wav=wav, toff=toff, nsmooth=nsmooth)

f1_inv_minus_f, f1_inv_plus_f, p0_minus_f, g_inv_minus_f, g_inv_plus_f = \
    MarchenkoWM.apply_onepoint(directVS, G0=G0sub.T, rtm=True, greens=True, 
                               dottest=True, **dict(iter_lim=n_iter, show=True))
g_inv_tot_f = g_inv_minus_f + g_inv_plus_f

Dot test passed, v^T(Opu)=-833.527331 - u^T(Op^Tv)=-833.527331
Dot test passed, v^T(Opu)=269.919973 - u^T(Op^Tv)=269.919973
 
LSQR            Least-squares solution of  Ax = b
The matrix A has   322998 rows  and   322998 cols
damp = 0.00000000000000e+00   calc_var =        0
atol = 1.00e-08                 conlim = 1.00e+08
btol = 1.00e-08               iter_lim =       15
 
   Itn      x[0]       r1norm     r2norm   Compatible    LS      Norm A   Cond A
     0  0.00000e+00   2.983e+07  2.983e+07    1.0e+00  3.5e-08
     1  0.00000e+00   1.311e+07  1.311e+07    4.4e-01  9.2e-01   1.1e+00  1.0e+00
     2  0.00000e+00   7.406e+06  7.406e+06    2.5e-01  3.9e-01   1.8e+00  2.2e+00
     3  0.00000e+00   5.479e+06  5.479e+06    1.8e-01  3.3e-01   2.1e+00  3.4e+00
     4  0.00000e+00   3.659e+06  3.659e+06    1.2e-01  3.3e-01   2.5e+00  5.2e+00
     5  0.00000e+00   2.780e+06  2.780e+06    9.3e-02  2.6e-01   2.9e+00  6.8e+00
     6  0.00000e+00   2.244e+06  2.244e+06    7.5e-02  2.3e-01   3.2e+00  8.5e+00
     7  0.00000e+00   1.498e+06  1.498e+06    5.0e-02  2.5e-01   3.6e+00  1.1e+01
     8  0.00000e+00   1.105e+06  1.105e+06    3.7e-02  1.9e-01   3.9e+00  1.3e+01
     9  0.00000e+00   8.988e+05  8.988e+05    3.0e-02  1.6e-01   4.2e+00  1.4e+01
    10  0.00000e+00   6.714e+05  6.714e+05    2.3e-02  1.7e-01   4.4e+00  1.6e+01
    11  0.00000e+00   5.040e+05  5.040e+05    1.7e-02  1.5e-01   4.6e+00  1.8e+01
    12  0.00000e+00   4.189e+05  4.189e+05    1.4e-02  1.4e-01   4.8e+00  2.0e+01
    13  0.00000e+00   3.357e+05  3.357e+05    1.1e-02  1.4e-01   5.0e+00  2.2e+01
    14  0.00000e+00   2.718e+05  2.718e+05    9.1e-03  1.1e-01   5.3e+00  2.4e+01
    15  0.00000e+00   2.424e+05  2.424e+05    8.1e-03  7.9e-02   5.5e+00  2.6e+01
 
LSQR finished
The iteration limit has been reached                      
 
istop =       7   r1norm = 2.4e+05   anorm = 5.5e+00   arnorm = 1.0e+05
itn   =      15   r2norm = 2.4e+05   acond = 2.6e+01   xnorm  = 3.5e+07
In [22]:
%timeit -r 3 -n 3 MarchenkoWM.apply_onepoint(directVS, G0=G0sub.T, **dict(iter_lim=1))

1.19 s ± 58.4 ms per loop (mean ± std. dev. of 3 runs, 3 loops each)
In [23]:
fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_minus_f.T, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].plot(r[0], directVS,'b'),axs[0].plot(r[0], -directVS,'b')
axs[0].axis('tight')
axs[0].set_ylim(1, -1);
axs[1].imshow(f1_inv_minus_f.T, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$f^-$'), axs[0].set_xlabel(r'$x_R$')
axs[1].plot(r[0], directVS,'b'),axs[1].plot(r[0], -directVS,'b')
axs[1].axis('tight')
axs[1].set_ylim(1, -1);
axs[2].imshow(f1_inv_plus_f.T, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$f^+$'), axs[0].set_xlabel(r'$x_R$')
axs[2].plot(r[0], directVS,'b'),axs[2].plot(r[0], -directVS,'b')
axs[2].axis('tight')
axs[2].set_ylim(1, -1);

fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_minus_f.T, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(2, 0);
axs[1].imshow(g_inv_minus_f.T, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$g^-$'), axs[0].set_xlabel(r'$x_R$')
axs[1].axis('tight')
axs[1].set_ylim(2, 0);
axs[2].imshow(g_inv_plus_f.T, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$g^+$'), axs[0].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(2, 0);

fig = plt.figure(figsize=(15,9))
ax1 = plt.subplot2grid((1, 5), (0, 0), colspan=2)
ax2 = plt.subplot2grid((1, 5), (0, 2), colspan=2)
ax3 = plt.subplot2grid((1, 5), (0, 4))

ax1.imshow(Gsub, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], t[0]))
ax1.set_title(r'$G_{true}$'), ax1.set_xlabel(r'$x_R$'), ax1.set_ylabel(r'$t$')
ax1.axis('tight')
ax1.set_ylim(1.2, 0)
ax2.imshow(g_inv_tot_f.T, cmap='gray', vmin=-5e5, vmax=5e5, 
              extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
ax2.set_title(r'$G_{est}$'), ax2.set_xlabel(r'$x_R$'), ax2.set_ylabel(r'$t$')
ax2.axis('tight')
ax2.set_ylim(1.2, 0)

ax3.plot(Gsub[:, nr//2]/Gsub.max(), t, 'r', lw=5)
ax3.plot(g_inv_tot_f[nr//2, nt-1:]/g_inv_tot_f.max(), t, 'k', lw=3)
ax3.set_ylim(1.2, 0);

We can also provide the reflection response in the original time domain

In [24]:
MarchenkoWM = Marchenko(R, dt=dt, dr=dr, nfmax=nfmax, wav=wav,
                        toff=toff, nsmooth=nsmooth)

f1_inv_minus_f, f1_inv_plus_f, p0_minus_f, g_inv_minus_f, g_inv_plus_f = \
    MarchenkoWM.apply_onepoint(directVS, G0=G0sub.T, rtm=True, greens=True, 
                               dottest=True, **dict(iter_lim=n_iter, show=True))
g_inv_tot_f = g_inv_minus_f + g_inv_plus_f

Dot test passed, v^T(Opu)=-583.620241 - u^T(Op^Tv)=-583.620241
Dot test passed, v^T(Opu)=20.160749 - u^T(Op^Tv)=20.160749
 
LSQR            Least-squares solution of  Ax = b
The matrix A has   322998 rows  and   322998 cols
damp = 0.00000000000000e+00   calc_var =        0
atol = 1.00e-08                 conlim = 1.00e+08
btol = 1.00e-08               iter_lim =       15
 
   Itn      x[0]       r1norm     r2norm   Compatible    LS      Norm A   Cond A
     0  0.00000e+00   2.983e+07  2.983e+07    1.0e+00  3.5e-08
     1  0.00000e+00   1.311e+07  1.311e+07    4.4e-01  9.2e-01   1.1e+00  1.0e+00
     2  0.00000e+00   7.406e+06  7.406e+06    2.5e-01  3.9e-01   1.8e+00  2.2e+00
     3  0.00000e+00   5.479e+06  5.479e+06    1.8e-01  3.3e-01   2.1e+00  3.4e+00
     4  0.00000e+00   3.659e+06  3.659e+06    1.2e-01  3.3e-01   2.5e+00  5.2e+00
     5  0.00000e+00   2.780e+06  2.780e+06    9.3e-02  2.6e-01   2.9e+00  6.8e+00
     6  0.00000e+00   2.244e+06  2.244e+06    7.5e-02  2.3e-01   3.2e+00  8.5e+00
     7  0.00000e+00   1.498e+06  1.498e+06    5.0e-02  2.5e-01   3.6e+00  1.1e+01
     8  0.00000e+00   1.105e+06  1.105e+06    3.7e-02  1.9e-01   3.9e+00  1.3e+01
     9  0.00000e+00   8.988e+05  8.988e+05    3.0e-02  1.6e-01   4.2e+00  1.4e+01
    10  0.00000e+00   6.714e+05  6.714e+05    2.3e-02  1.7e-01   4.4e+00  1.6e+01
    11  0.00000e+00   5.040e+05  5.040e+05    1.7e-02  1.5e-01   4.6e+00  1.8e+01
    12  0.00000e+00   4.189e+05  4.189e+05    1.4e-02  1.4e-01   4.8e+00  2.0e+01
    13  0.00000e+00   3.357e+05  3.357e+05    1.1e-02  1.4e-01   5.0e+00  2.2e+01
    14  0.00000e+00   2.718e+05  2.718e+05    9.1e-03  1.1e-01   5.3e+00  2.4e+01
    15  0.00000e+00   2.424e+05  2.424e+05    8.1e-03  7.9e-02   5.5e+00  2.6e+01
 
LSQR finished
The iteration limit has been reached                      
 
istop =       7   r1norm = 2.4e+05   anorm = 5.5e+00   arnorm = 1.0e+05
itn   =      15   r2norm = 2.4e+05   acond = 2.6e+01   xnorm  = 3.5e+07
In [25]:
fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_minus_f.T, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].plot(r[0], directVS,'b'),axs[0].plot(r[0], -directVS,'b')
axs[0].axis('tight')
axs[0].set_ylim(1, -1);
axs[1].imshow(f1_inv_minus_f.T, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$f^-$'), axs[0].set_xlabel(r'$x_R$')
axs[1].plot(r[0], directVS,'b'),axs[1].plot(r[0], -directVS,'b')
axs[1].axis('tight')
axs[1].set_ylim(1, -1);
axs[2].imshow(f1_inv_plus_f.T, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$f^+$'), axs[0].set_xlabel(r'$x_R$')
axs[2].plot(r[0], directVS,'b'),axs[2].plot(r[0], -directVS,'b')
axs[2].axis('tight')
axs[2].set_ylim(1, -1);

fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_minus_f.T, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(2, 0);
axs[1].imshow(g_inv_minus_f.T, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$g^-$'), axs[0].set_xlabel(r'$x_R$')
axs[1].axis('tight')
axs[1].set_ylim(2, 0);
axs[2].imshow(g_inv_plus_f.T, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$g^+$'), axs[0].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(2, 0);

fig = plt.figure(figsize=(15,9))
ax1 = plt.subplot2grid((1, 5), (0, 0), colspan=2)
ax2 = plt.subplot2grid((1, 5), (0, 2), colspan=2)
ax3 = plt.subplot2grid((1, 5), (0, 4))

ax1.imshow(Gsub, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], t[0]))
ax1.set_title(r'$G_{true}$'), ax1.set_xlabel(r'$x_R$'), ax1.set_ylabel(r'$t$')
ax1.axis('tight')
ax1.set_ylim(1.2, 0)
ax2.imshow(g_inv_tot_f.T, cmap='gray', vmin=-5e5, vmax=5e5, 
              extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
ax2.set_title(r'$G_{est}$'), ax2.set_xlabel(r'$x_R$'), ax2.set_ylabel(r'$t$')
ax2.axis('tight')
ax2.set_ylim(1.2, 0)

ax3.plot(Gsub[:, nr//2]/Gsub.max(), t, 'r', lw=5)
ax3.plot(g_inv_tot_f[nr//2, nt-1:]/g_inv_tot_f.max(), t, 'k', lw=3)
ax3.set_ylim(1.2, 0);

Similarly, but we can now use the flag fast=True to allow fast calculation of MDC (see documentation for more details). This approach is however more demanding in terms of memory and should not be used for very large datasets.

In [26]:
MarchenkoWM = Marchenko(R, dt=dt, dr=dr, nfmax=nfmax, wav=wav,
                        toff=toff, nsmooth=nsmooth)

f1_inv_minus_f, f1_inv_plus_f, p0_minus_f, g_inv_minus_f, g_inv_plus_f = \
    MarchenkoWM.apply_onepoint(directVS, G0=G0sub.T, rtm=True, greens=True, 
                               dottest=True, fast=True, **dict(iter_lim=n_iter, show=True))
g_inv_tot_f = g_inv_minus_f + g_inv_plus_f

Dot test passed, v^T(Opu)=1029.272052 - u^T(Op^Tv)=1029.272052
Dot test passed, v^T(Opu)=23.623077 - u^T(Op^Tv)=23.623077
 
LSQR            Least-squares solution of  Ax = b
The matrix A has   322998 rows  and   322998 cols
damp = 0.00000000000000e+00   calc_var =        0
atol = 1.00e-08                 conlim = 1.00e+08
btol = 1.00e-08               iter_lim =       15
 
   Itn      x[0]       r1norm     r2norm   Compatible    LS      Norm A   Cond A
     0  0.00000e+00   2.983e+07  2.983e+07    1.0e+00  3.5e-08
     1  0.00000e+00   1.311e+07  1.311e+07    4.4e-01  9.2e-01   1.1e+00  1.0e+00
     2  0.00000e+00   7.406e+06  7.406e+06    2.5e-01  3.9e-01   1.8e+00  2.2e+00
     3  0.00000e+00   5.479e+06  5.479e+06    1.8e-01  3.3e-01   2.1e+00  3.4e+00
     4  0.00000e+00   3.659e+06  3.659e+06    1.2e-01  3.3e-01   2.5e+00  5.2e+00
     5  0.00000e+00   2.780e+06  2.780e+06    9.3e-02  2.6e-01   2.9e+00  6.8e+00
     6  0.00000e+00   2.244e+06  2.244e+06    7.5e-02  2.3e-01   3.2e+00  8.5e+00
     7  0.00000e+00   1.498e+06  1.498e+06    5.0e-02  2.5e-01   3.6e+00  1.1e+01
     8  0.00000e+00   1.105e+06  1.105e+06    3.7e-02  1.9e-01   3.9e+00  1.3e+01
     9  0.00000e+00   8.988e+05  8.988e+05    3.0e-02  1.6e-01   4.2e+00  1.4e+01
    10  0.00000e+00   6.714e+05  6.714e+05    2.3e-02  1.7e-01   4.4e+00  1.6e+01
    11  0.00000e+00   5.040e+05  5.040e+05    1.7e-02  1.5e-01   4.6e+00  1.8e+01
    12  0.00000e+00   4.189e+05  4.189e+05    1.4e-02  1.4e-01   4.8e+00  2.0e+01
    13  0.00000e+00   3.357e+05  3.357e+05    1.1e-02  1.4e-01   5.0e+00  2.2e+01
    14  0.00000e+00   2.718e+05  2.718e+05    9.1e-03  1.1e-01   5.3e+00  2.4e+01
    15  0.00000e+00   2.424e+05  2.424e+05    8.1e-03  7.9e-02   5.5e+00  2.6e+01
 
LSQR finished
The iteration limit has been reached                      
 
istop =       7   r1norm = 2.4e+05   anorm = 5.5e+00   arnorm = 1.0e+05
itn   =      15   r2norm = 2.4e+05   acond = 2.6e+01   xnorm  = 3.5e+07
In [27]:
fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_minus_f.T, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].plot(r[0], directVS,'b'),axs[0].plot(r[0], -directVS,'b')
axs[0].axis('tight')
axs[0].set_ylim(1, -1);
axs[1].imshow(f1_inv_minus_f.T, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$f^-$'), axs[0].set_xlabel(r'$x_R$')
axs[1].plot(r[0], directVS,'b'),axs[1].plot(r[0], -directVS,'b')
axs[1].axis('tight')
axs[1].set_ylim(1, -1);
axs[2].imshow(f1_inv_plus_f.T, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$f^+$'), axs[0].set_xlabel(r'$x_R$')
axs[2].plot(r[0], directVS,'b'),axs[2].plot(r[0], -directVS,'b')
axs[2].axis('tight')
axs[2].set_ylim(1, -1);

fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_minus_f.T, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(2, 0);
axs[1].imshow(g_inv_minus_f.T, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$g^-$'), axs[0].set_xlabel(r'$x_R$')
axs[1].axis('tight')
axs[1].set_ylim(2, 0);
axs[2].imshow(g_inv_plus_f.T, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$g^+$'), axs[0].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(2, 0);

fig = plt.figure(figsize=(15,9))
ax1 = plt.subplot2grid((1, 5), (0, 0), colspan=2)
ax2 = plt.subplot2grid((1, 5), (0, 2), colspan=2)
ax3 = plt.subplot2grid((1, 5), (0, 4))

ax1.imshow(Gsub, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], t[0]))
ax1.set_title(r'$G_{true}$'), ax1.set_xlabel(r'$x_R$'), ax1.set_ylabel(r'$t$')
ax1.axis('tight')
ax1.set_ylim(1.2, 0)
ax2.imshow(g_inv_tot_f.T, cmap='gray', vmin=-5e5, vmax=5e5, 
              extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
ax2.set_title(r'$G_{est}$'), ax2.set_xlabel(r'$x_R$'), ax2.set_ylabel(r'$t$')
ax2.axis('tight')
ax2.set_ylim(1.2, 0)

ax3.plot(Gsub[:, nr//2]/Gsub.max(), t, 'r', lw=5)
ax3.plot(g_inv_tot_f[nr//2, nt-1:]/g_inv_tot_f.max(), t, 'k', lw=3)
ax3.set_ylim(1.2, 0);

Finally we show that is also possible to use an analytical direct wave

In [28]:
MarchenkoWM = Marchenko(R, dt=dt, dr=dr, nfmax=nfmax, wav=wav,
                        toff=toff, nsmooth=nsmooth)

f1_inv_minus_f, f1_inv_plus_f, p0_minus_f, g_inv_minus_f, g_inv_plus_f = \
    MarchenkoWM.apply_onepoint(directVS, nfft=2**11, rtm=True, greens=True, 
                               dottest=True, **dict(iter_lim=n_iter, show=True))
g_inv_tot_f = g_inv_minus_f + g_inv_plus_f

Dot test passed, v^T(Opu)=-417.756482 - u^T(Op^Tv)=-417.756482
Dot test passed, v^T(Opu)=-88.329503 - u^T(Op^Tv)=-88.329503
 
LSQR            Least-squares solution of  Ax = b
The matrix A has   322998 rows  and   322998 cols
damp = 0.00000000000000e+00   calc_var =        0
atol = 1.00e-08                 conlim = 1.00e+08
btol = 1.00e-08               iter_lim =       15
 
   Itn      x[0]       r1norm     r2norm   Compatible    LS      Norm A   Cond A
     0  0.00000e+00   3.945e+01  3.945e+01    1.0e+00  2.6e-02
     1  0.00000e+00   1.734e+01  1.734e+01    4.4e-01  9.2e-01   1.1e+00  1.0e+00
     2  0.00000e+00   9.792e+00  9.792e+00    2.5e-01  3.9e-01   1.8e+00  2.2e+00
     3  0.00000e+00   7.242e+00  7.242e+00    1.8e-01  3.3e-01   2.1e+00  3.4e+00
     4  0.00000e+00   4.838e+00  4.838e+00    1.2e-01  3.3e-01   2.5e+00  5.2e+00
     5  0.00000e+00   3.674e+00  3.674e+00    9.3e-02  2.6e-01   2.9e+00  6.8e+00
     6  0.00000e+00   2.963e+00  2.963e+00    7.5e-02  2.3e-01   3.2e+00  8.5e+00
     7  0.00000e+00   1.980e+00  1.980e+00    5.0e-02  2.5e-01   3.6e+00  1.1e+01
     8  0.00000e+00   1.460e+00  1.460e+00    3.7e-02  1.9e-01   3.9e+00  1.3e+01
     9  0.00000e+00   1.187e+00  1.187e+00    3.0e-02  1.6e-01   4.2e+00  1.4e+01
    10  0.00000e+00   8.854e-01  8.854e-01    2.2e-02  1.7e-01   4.4e+00  1.6e+01
    11  0.00000e+00   6.633e-01  6.633e-01    1.7e-02  1.5e-01   4.6e+00  1.8e+01
    12  0.00000e+00   5.497e-01  5.497e-01    1.4e-02  1.4e-01   4.8e+00  2.0e+01
    13  0.00000e+00   4.388e-01  4.388e-01    1.1e-02  1.4e-01   5.0e+00  2.2e+01
    14  0.00000e+00   3.534e-01  3.534e-01    9.0e-03  1.1e-01   5.3e+00  2.4e+01
    15  0.00000e+00   3.137e-01  3.137e-01    8.0e-03  8.0e-02   5.5e+00  2.6e+01
 
LSQR finished
The iteration limit has been reached                      
 
istop =       7   r1norm = 3.1e-01   anorm = 5.5e+00   arnorm = 1.4e-01
itn   =      15   r2norm = 3.1e-01   acond = 2.6e+01   xnorm  = 4.6e+01
In [29]:
fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_minus_f.T, cmap='gray', vmin=-5e-1, vmax=5e-1, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].plot(r[0], directVS,'b'),axs[0].plot(r[0], -directVS,'b')
axs[0].axis('tight')
axs[0].set_ylim(1, -1);
axs[1].imshow(f1_inv_minus_f.T, cmap='gray', vmin=-5e-1, vmax=5e-1, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$f^-$'), axs[0].set_xlabel(r'$x_R$')
axs[1].plot(r[0], directVS,'b'),axs[1].plot(r[0], -directVS,'b')
axs[1].axis('tight')
axs[1].set_ylim(1, -1);
axs[2].imshow(f1_inv_plus_f.T, cmap='gray', vmin=-5e-1, vmax=5e-1, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$f^+$'), axs[0].set_xlabel(r'$x_R$')
axs[2].plot(r[0], directVS,'b'),axs[2].plot(r[0], -directVS,'b')
axs[2].axis('tight')
axs[2].set_ylim(1, -1);

fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_minus_f.T, cmap='gray', vmin=-5e-1, vmax=5e-1, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(2, 0);
axs[1].imshow(g_inv_minus_f.T, cmap='gray', vmin=-5e-1, vmax=5e-1, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$g^-$'), axs[0].set_xlabel(r'$x_R$')
axs[1].axis('tight')
axs[1].set_ylim(2, 0);
axs[2].imshow(g_inv_plus_f.T, cmap='gray', vmin=-5e-1, vmax=5e-1, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$g^+$'), axs[0].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(2, 0);

fig = plt.figure(figsize=(15,9))
ax1 = plt.subplot2grid((1, 5), (0, 0), colspan=2)
ax2 = plt.subplot2grid((1, 5), (0, 2), colspan=2)
ax3 = plt.subplot2grid((1, 5), (0, 4))

ax1.imshow(Gsub, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], t[0]))
ax1.set_title(r'$G_{true}$'), ax1.set_xlabel(r'$x_R$'), ax1.set_ylabel(r'$t$')
ax1.axis('tight')
ax1.set_ylim(1.2, 0)
ax2.imshow(g_inv_tot_f.T, cmap='gray', vmin=-5e-1, vmax=5e-1,
              extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
ax2.set_title(r'$G_{est}$'), ax2.set_xlabel(r'$x_R$'), ax2.set_ylabel(r'$t$')
ax2.axis('tight')
ax2.set_ylim(1.2, 0)

ax3.plot(Gsub[:, nr//2]/Gsub.max(), t, 'r', lw=5)
ax3.plot(g_inv_tot_f[nr//2, nt-1:]/g_inv_tot_f.max(), t, 'k', lw=3)
ax3.set_ylim(1.2, 0);

Examples of Marchenko class usage for multiple points

Finally, lets consider an example with multiple virtual points

In [33]:
nvs = 11
vs = [np.arange(11)*100 + 1000, 
      np.ones(11)*1060]

plt.figure(figsize=(18,9))
plt.imshow(rho, cmap='gray', extent = (x[0], x[-1], z[-1], z[0]))
plt.scatter(s[0, 5::10], s[1, 5::10], marker='*', s=150, c='r', edgecolors='k')
plt.scatter(r[0, ::10],  r[1, ::10], marker='v', s=150, c='b', edgecolors='k')
plt.scatter(vs[0], vs[1], marker='.', s=250, c='m', edgecolors='k')
plt.axis('tight')
plt.xlabel('x [m]'),plt.ylabel('y [m]'),plt.title('Model and Geometry')
plt.xlim(x[0], x[-1]);

Compute direct arrival for all pairs of receivers and virtual sources

In [34]:
# direct arrival window - traveltime
directVS = np.sqrt((vs[0]-r[0][:, np.newaxis])**2+(vs[1]-r[1][:, np.newaxis])**2)/vel
directVS_off = directVS - toff

plt.figure()
im = plt.imshow(directVS, cmap='gist_rainbow')
plt.axis('tight')
plt.xlabel('#VS'),plt.ylabel('#Rec'),plt.title('Direct arrival')
plt.colorbar(im);

Create the window and analytical direct arrival

In [35]:
idirectVS_off = np.round(directVS_off/dt).astype(np.int)
w = np.zeros((nr, nvs, nt))
for ir in range(nr):
    for ivs in range(nvs):
        w[ir, ivs, :idirectVS_off[ir, ivs]]=1            
w = np.concatenate((np.flip(w, axis=-1), w[:,:, 1:]), axis=-1)

if nsmooth>0:
    smooth=np.ones(nsmooth)/nsmooth
    w  = filtfilt(smooth, 1, w)    

fig, axs = plt.subplots(1, 2,  sharey=True, figsize=(15, 9))
im = axs[0].imshow(w[nr//2].T, cmap='gray', extent=(vs[0][0], vs[0][-1], t[-1], -t[-1]))
axs[0].plot(vs[0], directVS[nr//2],'b'),axs[0].plot(vs[0], -directVS[nr//2],'b')
axs[0].set_title('Window (fixed r)'), axs[0].set_xlabel(r'$x_{VS}$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
fig.colorbar(im, ax=axs[0])
im = axs[1].imshow(w[:, nvs//2].T, cmap='gray', 
                   extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].plot(r[0], directVS[:, nvs//2],'b'),axs[1].plot(r[0], -directVS[:, nvs//2],'b')
axs[1].set_title('Window (fixed vs)'), axs[1].set_xlabel(r'$x_R$'), axs[1].set_ylabel(r'$t$')
axs[1].axis('tight')
fig.colorbar(im, ax=axs[1]);
In [36]:
G0sub_ana = np.zeros((nr, nvs, nt))
for ivs in range(nvs):
    G0sub_ana[:, ivs] = directwave(wav, directVS[:,ivs], nt, dt, nfft=int(2**(np.ceil(np.log2(nt))))).T

fig, axs = plt.subplots(1, 2,  sharey=True, figsize=(15, 9))
im = axs[0].imshow(G0sub_ana[nr//2].T, cmap='gray', extent=(vs[0][0], vs[0][-1], t[-1], t[0]))
axs[0].plot(vs[0], directVS[nr//2],'b')
axs[0].set_title('G0 (fixed r)'), axs[0].set_xlabel(r'$x_{VS}$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
fig.colorbar(im, ax=axs[0])
im = axs[1].imshow(G0sub_ana[:,nvs//2].T, cmap='gray', extent=(r[0,0], r[0,-1], t[-1], t[0]))
axs[1].plot(r[0], directVS[:, nvs//2],'b')
axs[1].set_title('G0 (fixed vs)'), axs[1].set_xlabel(r'$x_R$'), axs[1].set_ylabel(r'$t$')
axs[1].axis('tight')
fig.colorbar(im, ax=axs[1]);

Let's now create and apply the Marchenko operator as done above for single virtual point

In [37]:
# Add negative time to operators
Rtwosided = np.concatenate((np.zeros((nr, ns, nt-1)), R), axis=-1)
R1twosided = np.concatenate((np.swapaxes(np.fliplr(np.swapaxes(R, 1, 2)), 1, 2), 
                            np.zeros((nr, ns, nt-1))), axis=-1)

Rtwosided_fft = np.fft.rfft(Rtwosided, 2*nt-1, axis=-1)/np.sqrt(2*nt-1)
Rtwosided_fft = Rtwosided_fft[...,:nfmax]
R1twosided_fft = np.fft.rfft(R1twosided, 2*nt-1, axis=-1)/np.sqrt(2*nt-1)
R1twosided_fft = R1twosided_fft[...,:nfmax]

# Operators
Rop = MDC(Rtwosided_fft, nt=2*nt-1, nv=nvs, dt=dt, dr=ds, twosided=True, dtype='complex64')
dottest(Rop, (2*nt-1)*nr*nvs, (2*nt-1)*ns*nvs, verb=True)

R1op = MDC(R1twosided_fft, nt=2*nt-1, nv=nvs, dt=dt, dr=ds, twosided=True, dtype='complex64')
dottest(R1op, (2*nt-1)*nr*nvs, (2*nt-1)*ns*nvs, verb=True)

# Input focusing function
fd_plus =  np.concatenate((np.flip(G0sub_ana, axis=-1), 
                           np.zeros((nr, nvs, nt-1))), axis=-1)

Dot test passed, v^T(Opu)=-301.872405 - u^T(Op^Tv)=-301.872405
Dot test passed, v^T(Opu)=25.652942 - u^T(Op^Tv)=25.652942
In [38]:
Wop = Diagonal(w.flatten())
Iop = Identity(nr*nvs*(2*nt-1))

Mop = VStack([HStack([Iop, -1*Wop*Rop]),
              HStack([-1*Wop*R1op, Iop])])*BlockDiag([Wop, Wop])
Gop = VStack([HStack([Iop, -1*Rop]),
              HStack([-1*R1op, Iop])])

dottest(Gop, 2*nr*nvs*(2*nt-1), 2*nr*nvs*(2*nt-1), verb=True)
dottest(Mop, 2*nr*nvs*(2*nt-1), 2*nr*nvs*(2*nt-1), verb=True);

p0_minus = Rop * fd_plus.flatten()
p0_minus = p0_minus.reshape(nr, nvs, (2*nt-1))

p0_minus = Rop * fd_plus.flatten()
p0_minus = p0_minus.reshape(nr, nvs, (2*nt-1))

d = Wop*Rop*fd_plus.flatten()
d = np.concatenate((d.reshape(nr, nvs, 2*nt-1), np.zeros((nr, nvs, 2*nt-1))))

f1_adj = Mop.H*d.flatten()
f1_inv = lsqr(Mop, d.flatten(), iter_lim=n_iter, show=True)[0]

f1_adj = f1_adj.reshape(2*nr, nvs, (2*nt-1))
f1_inv = f1_inv.reshape(2*nr, nvs, (2*nt-1))

Dot test passed, v^T(Opu)=21.979308 - u^T(Op^Tv)=21.979308
Dot test passed, v^T(Opu)=239.482065 - u^T(Op^Tv)=239.482065
 
LSQR            Least-squares solution of  Ax = b
The matrix A has 3.55298e+06 rows  and 3.55298e+06 cols
damp = 0.00000000000000e+00   calc_var =        0
atol = 1.00e-08                 conlim = 1.00e+08
btol = 1.00e-08               iter_lim =       15
 
   Itn      x[0]       r1norm     r2norm   Compatible    LS      Norm A   Cond A
     0  0.00000e+00   1.306e+02  1.306e+02    1.0e+00  7.9e-03
     1  0.00000e+00   5.675e+01  5.675e+01    4.3e-01  9.4e-01   1.1e+00  1.0e+00
     2  0.00000e+00   3.199e+01  3.199e+01    2.4e-01  3.9e-01   1.8e+00  2.2e+00
     3  0.00000e+00   2.376e+01  2.376e+01    1.8e-01  3.3e-01   2.1e+00  3.4e+00
     4  0.00000e+00   1.613e+01  1.613e+01    1.2e-01  3.2e-01   2.5e+00  5.1e+00
     5  0.00000e+00   1.244e+01  1.244e+01    9.5e-02  2.6e-01   2.9e+00  6.8e+00
     6  0.00000e+00   1.003e+01  1.003e+01    7.7e-02  2.2e-01   3.3e+00  8.5e+00
     7  0.00000e+00   6.938e+00  6.938e+00    5.3e-02  2.4e-01   3.6e+00  1.1e+01
     8  0.00000e+00   5.215e+00  5.215e+00    4.0e-02  1.9e-01   3.9e+00  1.3e+01
     9  0.00000e+00   4.198e+00  4.198e+00    3.2e-02  1.6e-01   4.2e+00  1.4e+01
    10  0.00000e+00   3.162e+00  3.162e+00    2.4e-02  1.8e-01   4.4e+00  1.6e+01
    11  0.00000e+00   2.392e+00  2.392e+00    1.8e-02  1.6e-01   4.6e+00  1.8e+01
    12  0.00000e+00   1.957e+00  1.957e+00    1.5e-02  1.3e-01   4.9e+00  2.0e+01
    13  0.00000e+00   1.542e+00  1.542e+00    1.2e-02  1.5e-01   5.0e+00  2.2e+01
    14  0.00000e+00   1.241e+00  1.241e+00    9.5e-03  1.2e-01   5.3e+00  2.4e+01
    15  0.00000e+00   1.093e+00  1.093e+00    8.4e-03  7.9e-02   5.5e+00  2.6e+01
 
LSQR finished
The iteration limit has been reached                      
 
istop =       7   r1norm = 1.1e+00   anorm = 5.5e+00   arnorm = 4.8e-01
itn   =      15   r2norm = 1.1e+00   acond = 2.6e+01   xnorm  = 1.5e+02
In [39]:
f1_adj_tot = f1_adj + np.concatenate((np.zeros((nr, nvs, 2*nt-1)),
                                      fd_plus))

f1_inv_tot = f1_inv + np.concatenate((np.zeros((nr, nvs, 2*nt-1)),
                                      fd_plus))

g_adj = Gop*f1_adj_tot.flatten()
g_inv = Gop*f1_inv_tot.flatten()

g_adj = g_adj.reshape(2*nr, nvs, (2*nt-1))
g_inv = g_inv.reshape(2*nr, nvs, (2*nt-1))

f1_adj_minus, f1_adj_plus =  f1_adj_tot[:nr], f1_adj_tot[nr:]
f1_inv_minus, f1_inv_plus =  f1_inv_tot[:nr], f1_inv_tot[nr:]

g_adj_minus, g_adj_plus =  -g_adj[:nr], np.flip(g_adj[nr:], axis=-1)
g_inv_minus, g_inv_plus =  -g_inv[:nr], np.flip(g_inv[nr:], axis=-1)

g_inv_tot = g_inv_minus + g_inv_plus
In [40]:
fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_minus[:, 2].T, cmap='gray', 
              vmin=-5e-1, vmax=5e-1, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].plot(r[0], directVS[:, 2],'b'),axs[0].plot(r[0], -directVS[:, 2],'b')
axs[0].axis('tight')
axs[0].set_ylim(1, -1);
axs[1].imshow(f1_inv_minus[:, 2].T, cmap='gray', vmin=-5e-1, vmax=5e-1, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$f^-$'), axs[0].set_xlabel(r'$x_R$')
axs[1].plot(r[0], directVS[:, 2],'b'),axs[1].plot(r[0], -directVS[:, 2],'b')
axs[1].axis('tight')
axs[1].set_ylim(1, -1);
axs[2].imshow(f1_inv_plus[:, 2].T, cmap='gray', vmin=-5e-1, vmax=5e-1, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$f^+$'), axs[0].set_xlabel(r'$x_R$')
axs[2].plot(r[0], directVS[:, 2],'b'),axs[2].plot(r[0], -directVS[:, 2],'b')
axs[2].axis('tight')
axs[2].set_ylim(1, -1);

fig, axs = plt.subplots(1, 3, sharey=True, figsize=(16, 9))
axs[0].imshow(p0_minus[:, 2].T, cmap='gray', vmin=-5e-1, vmax=5e-1, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(2, 0);
axs[1].imshow(g_inv_minus[:, 2].T, cmap='gray', vmin=-5e-1, vmax=5e-1, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[1].set_title(r'$g^-$'), axs[0].set_xlabel(r'$x_R$')
axs[1].axis('tight')
axs[1].set_ylim(2, 0);
axs[2].imshow(g_inv_plus[:, 2].T, cmap='gray', vmin=-5e-1, vmax=5e-1, extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
axs[2].set_title(r'$g^+$'), axs[0].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(2, 0);

fig = plt.figure(figsize=(15,9))
ax1 = plt.subplot2grid((1, 5), (0, 0), colspan=2)
ax2 = plt.subplot2grid((1, 5), (0, 2), colspan=2)
ax3 = plt.subplot2grid((1, 5), (0, 4))

ax1.imshow(Gsub, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], t[0]))
ax1.set_title(r'$G_{true}$'), ax1.set_xlabel(r'$x_R$'), ax1.set_ylabel(r'$t$')
ax1.axis('tight')
ax1.set_ylim(1.2, 0)
ax2.imshow(g_inv_tot[:, 2].T, cmap='gray', vmin=-5e-1, vmax=5e-1, 
              extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
ax2.set_title(r'$G_{est}$'), ax2.set_xlabel(r'$x_R$'), ax2.set_ylabel(r'$t$')
ax2.axis('tight')
ax2.set_ylim(1.2, 0)

ax3.plot(Gsub[:, nr//2]/Gsub.max(), t, 'r', lw=5)
ax3.plot(g_inv_tot[nr//2, 2, nt-1:]/g_inv_tot.max(), t, 'k', lw=3)
ax3.set_ylim(1.2, 0);

Finally we use the apply_multiplepoints function to compute the response to multiple virtual points in the subsurface

In [41]:
MarchenkoWM = Marchenko(R, dt=dt, dr=dr, nfmax=nfmax, wav=wav,
                        toff=toff, nsmooth=nsmooth)
In [42]:
f1_inv_minus_f, f1_inv_plus_f, p0_minus_f, g_inv_minus_f, g_inv_plus_f = \
    MarchenkoWM.apply_multiplepoints(directVS, nfft=2**11, rtm=True, greens=True, 
                               dottest=False, **dict(iter_lim=n_iter, show=True))
g_inv_tot_f = g_inv_minus_f + g_inv_plus_f

 
LSQR            Least-squares solution of  Ax = b
The matrix A has 3.55298e+06 rows  and 3.55298e+06 cols
damp = 0.00000000000000e+00   calc_var =        0
atol = 1.00e-08                 conlim = 1.00e+08
btol = 1.00e-08               iter_lim =       15
 
   Itn      x[0]       r1norm     r2norm   Compatible    LS      Norm A   Cond A
     0  0.00000e+00   1.306e+02  1.306e+02    1.0e+00  7.9e-03
     1  0.00000e+00   5.675e+01  5.675e+01    4.3e-01  9.4e-01   1.1e+00  1.0e+00
     2  0.00000e+00   3.199e+01  3.199e+01    2.4e-01  3.9e-01   1.8e+00  2.2e+00
     3  0.00000e+00   2.376e+01  2.376e+01    1.8e-01  3.3e-01   2.1e+00  3.4e+00
     4  0.00000e+00   1.613e+01  1.613e+01    1.2e-01  3.2e-01   2.5e+00  5.1e+00
     5  0.00000e+00   1.244e+01  1.244e+01    9.5e-02  2.6e-01   2.9e+00  6.8e+00
     6  0.00000e+00   1.003e+01  1.003e+01    7.7e-02  2.2e-01   3.3e+00  8.5e+00
     7  0.00000e+00   6.938e+00  6.938e+00    5.3e-02  2.4e-01   3.6e+00  1.1e+01
     8  0.00000e+00   5.215e+00  5.215e+00    4.0e-02  1.9e-01   3.9e+00  1.3e+01
     9  0.00000e+00   4.198e+00  4.198e+00    3.2e-02  1.6e-01   4.2e+00  1.4e+01
    10  0.00000e+00   3.162e+00  3.162e+00    2.4e-02  1.8e-01   4.4e+00  1.6e+01
    11  0.00000e+00   2.392e+00  2.392e+00    1.8e-02  1.6e-01   4.6e+00  1.8e+01
    12  0.00000e+00   1.957e+00  1.957e+00    1.5e-02  1.3e-01   4.9e+00  2.0e+01
    13  0.00000e+00   1.542e+00  1.542e+00    1.2e-02  1.5e-01   5.0e+00  2.2e+01
    14  0.00000e+00   1.241e+00  1.241e+00    9.5e-03  1.2e-01   5.3e+00  2.4e+01
    15  0.00000e+00   1.093e+00  1.093e+00    8.4e-03  7.9e-02   5.5e+00  2.6e+01
 
LSQR finished
The iteration limit has been reached                      
 
istop =       7   r1norm = 1.1e+00   anorm = 5.5e+00   arnorm = 4.8e-01
itn   =      15   r2norm = 1.1e+00   acond = 2.6e+01   xnorm  = 1.5e+02
In [43]:
fig, axs = plt.subplots(5, 1, figsize=(16, 22))
axs[0].imshow(np.swapaxes(p0_minus_f, 0, 1).reshape(nr*nvs, 2*nt-1).T, cmap='gray', 
              vmin=-5e-1, vmax=5e-1, extent=(0, nr*nvs, t[-1], -t[-1]))
axs[0].set_title(r'$p_0^-$'), axs[0].set_xlabel(r'$x_R$'), axs[0].set_ylabel(r'$t$')
axs[0].axis('tight')
axs[0].set_ylim(1, -1);
axs[1].imshow(np.swapaxes(f1_inv_minus_f, 0, 1).reshape(nr*nvs,2*nt-1).T, cmap='gray', 
              vmin=-5e-1, vmax=5e-1, extent=(0, nr*nvs, t[-1], -t[-1]))
axs[1].set_title(r'$f^-$'), axs[0].set_xlabel(r'$x_R$')
axs[1].axis('tight')
axs[1].set_ylim(1, -1);
axs[2].imshow(np.swapaxes(f1_inv_plus_f, 0, 1).reshape(nr*nvs,2*nt-1).T, cmap='gray', 
              vmin=-5e-1, vmax=5e-1, extent=(0, nr*nvs, t[-1], -t[-1]))
axs[2].set_title(r'$f^+$'), axs[0].set_xlabel(r'$x_R$')
axs[2].axis('tight')
axs[2].set_ylim(1, -1);
axs[3].imshow(np.swapaxes(g_inv_minus, 0, 1).reshape(nr*nvs,2*nt-1).T, cmap='gray', 
              vmin=-5e-1, vmax=5e-1, extent=(0, nr*nvs, t[-1], -t[-1]))
axs[3].set_title(r'$g^-$'), axs[0].set_xlabel(r'$x_R$')
axs[3].axis('tight')
axs[3].set_ylim(2, 0);
axs[4].imshow(np.swapaxes(g_inv_plus, 0, 1).reshape(nr*nvs,2*nt-1).T, cmap='gray', 
              vmin=-5e-1, vmax=5e-1, extent=(0, nr*nvs, t[-1], -t[-1]))
axs[4].set_title(r'$g^+$'), axs[0].set_xlabel(r'$x_R$')
axs[4].axis('tight')
axs[4].set_ylim(2, 0);
In [44]:
fig = plt.figure(figsize=(15,9))
ax1 = plt.subplot2grid((1, 5), (0, 0), colspan=2)
ax2 = plt.subplot2grid((1, 5), (0, 2), colspan=2)
ax3 = plt.subplot2grid((1, 5), (0, 4))

ax1.imshow(Gsub, cmap='gray', vmin=-5e5, vmax=5e5, extent=(r[0,0], r[0,-1], t[-1], t[0]))
ax1.set_title(r'$G_{true}$'), ax1.set_xlabel(r'$x_R$'), ax1.set_ylabel(r'$t$')
ax1.axis('tight')
ax1.set_ylim(1.2, 0)
ax2.imshow(g_inv_tot_f[:, 2].T, cmap='gray', vmin=-5e-1, vmax=5e-1, 
              extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
ax2.set_title(r'$G_{est}$'), ax2.set_xlabel(r'$x_R$'), ax2.set_ylabel(r'$t$')
ax2.axis('tight')
ax2.set_ylim(1.2, 0)

ax3.plot(Gsub[:, nr//2]/Gsub.max(), t, 'r', lw=5)
ax3.plot(g_inv_tot_f[nr//2, 2, nt-1:]/g_inv_tot_f.max(), t, 'k', lw=3)
ax3.set_ylim(1.2, 0);

Let's try now to do the same for a depth level $z_1=z_0 + dz$. We try to use the focusing functions at the level above as a starting guess

In [45]:
nvs = 11
dz = 50
vs = [np.arange(11)*100 + 1000, 
      np.ones(11)*1060 + dz]

# direct arrival window - traveltime
directVS = np.sqrt((vs[0]-r[0][:, np.newaxis])**2+(vs[1]-r[1][:, np.newaxis])**2)/vel
directVS_off = directVS - toff
idirectVS_off = np.round(directVS_off/dt).astype(np.int)
w = np.zeros((nr, nvs, nt))
for ir in range(nr):
    for ivs in range(nvs):
        w[ir, ivs, :idirectVS_off[ir, ivs]]=1            
w = np.concatenate((np.flip(w, axis=-1), w[:,:, 1:]), axis=-1)

if nsmooth>0:
    smooth=np.ones(nsmooth)/nsmooth
    w  = filtfilt(smooth, 1, w)    

G0sub_ana = np.zeros((nr, nvs, nt))
for ivs in range(nvs):
    G0sub_ana[:, ivs] = directwave(wav, directVS[:,ivs], nt, dt, nfft=int(2**(np.ceil(np.log2(nt))))).T

# Input focusing function
fd_plus_z1 =  np.concatenate((np.flip(G0sub_ana, axis=-1), 
                           np.zeros((nr, nvs, nt-1))), axis=-1)

Wop = Diagonal(w.flatten())
Iop = Identity(nr*nvs*(2*nt-1))

Mop = VStack([HStack([Iop, -1*Wop*Rop]),
              HStack([-1*Wop*R1op, Iop])])*BlockDiag([Wop, Wop])
Gop = VStack([HStack([Iop, -1*Rop]),
              HStack([-1*R1op, Iop])])

dottest(Gop, 2*nr*nvs*(2*nt-1), 2*nr*nvs*(2*nt-1), verb=True)
dottest(Mop, 2*nr*nvs*(2*nt-1), 2*nr*nvs*(2*nt-1), verb=True);

p0_minus = Rop * fd_plus_z1.flatten()
p0_minus = p0_minus.reshape(nr, nvs, (2*nt-1))

d = Wop*Rop*fd_plus_z1.flatten()
d = np.concatenate((d.reshape(nr, nvs, 2*nt-1), np.zeros((nr, nvs, 2*nt-1))))

# initial guess 
fstart = np.concatenate((f1_inv_minus_f, f1_inv_plus_f - fd_plus), axis=0)

# invert focusing
f1_inv_z1fast = RegularizedInversion(Mop, None, d.flatten(), x0=fstart.ravel(), 
                                     **dict(iter_lim=n_iter-2, show=True))
f1_inv_z1fast = f1_inv_z1fast.reshape(2*nr, nvs, (2*nt-1))
f1_inv_tot_z1fast = f1_inv_z1fast + np.concatenate((np.zeros((nr, nvs, 2*nt-1)),
                                                    fd_plus_z1))

g_inv_z1fast = Gop*f1_inv_tot_z1fast.flatten()
g_inv_z1fast = g_inv_z1fast.reshape(2*nr, nvs, (2*nt-1))
f1_inv_minus_z1fast, f1_inv_plus_z1fast =  f1_inv_tot_z1fast[:nr], f1_inv_tot_z1fast[nr:]

g_inv_minus_z1fast, g_inv_plus_z1fast = -g_inv_z1fast[:nr], np.flip(g_inv_z1fast[nr:], axis=-1)
g_inv_tot_z1fast = g_inv_minus_z1fast + g_inv_plus_z1fast

Dot test passed, v^T(Opu)=-1270.342318 - u^T(Op^Tv)=-1270.342318
Dot test passed, v^T(Opu)=-1155.882028 - u^T(Op^Tv)=-1155.882028
 
LSQR            Least-squares solution of  Ax = b
The matrix A has 3.55298e+06 rows  and 3.55298e+06 cols
damp = 0.00000000000000e+00   calc_var =        0
atol = 1.00e-08                 conlim = 1.00e+08
btol = 1.00e-08               iter_lim =       13
 
   Itn      x[0]       r1norm     r2norm   Compatible    LS      Norm A   Cond A
     0  0.00000e+00   2.299e+02  2.299e+02    1.0e+00  4.5e-03
     1  0.00000e+00   1.001e+02  1.001e+02    4.4e-01  9.4e-01   1.1e+00  1.0e+00
     2  0.00000e+00   5.638e+01  5.638e+01    2.5e-01  3.9e-01   1.8e+00  2.2e+00
     3  0.00000e+00   4.209e+01  4.209e+01    1.8e-01  3.2e-01   2.2e+00  3.4e+00
     4  0.00000e+00   2.869e+01  2.869e+01    1.2e-01  3.1e-01   2.5e+00  5.1e+00
     5  0.00000e+00   2.141e+01  2.141e+01    9.3e-02  2.6e-01   2.8e+00  6.8e+00
     6  0.00000e+00   1.745e+01  1.745e+01    7.6e-02  2.4e-01   3.2e+00  8.5e+00
     7  0.00000e+00   1.212e+01  1.212e+01    5.3e-02  2.3e-01   3.6e+00  1.1e+01
     8  0.00000e+00   9.197e+00  9.197e+00    4.0e-02  1.9e-01   3.8e+00  1.3e+01
     9  0.00000e+00   7.583e+00  7.583e+00    3.3e-02  1.6e-01   4.1e+00  1.4e+01
    10  0.00000e+00   5.700e+00  5.700e+00    2.5e-02  1.8e-01   4.4e+00  1.6e+01
    11  0.00000e+00   4.202e+00  4.202e+00    1.8e-02  1.6e-01   4.6e+00  1.8e+01
    12  0.00000e+00   3.383e+00  3.383e+00    1.5e-02  1.4e-01   4.9e+00  2.0e+01
    13  0.00000e+00   2.558e+00  2.558e+00    1.1e-02  1.6e-01   5.1e+00  2.2e+01
 
LSQR finished
The iteration limit has been reached                      
 
istop =       7   r1norm = 2.6e+00   anorm = 5.1e+00   arnorm = 2.1e+00
itn   =      13   r2norm = 2.6e+00   acond = 2.2e+01   xnorm  = 2.7e+02

Let's compare it without an initial guess

In [46]:
f1_inv_minus_z1, f1_inv_plus_z1, p0_minus_z1, g_inv_minus_z1, g_inv_plus_z1 = \
    MarchenkoWM.apply_multiplepoints(directVS, nfft=2**11, rtm=True, greens=True, 
                               dottest=False, **dict(iter_lim=n_iter-2, show=True))
g_inv_tot_z1 = g_inv_minus_z1 + g_inv_plus_z1

 
LSQR            Least-squares solution of  Ax = b
The matrix A has 3.55298e+06 rows  and 3.55298e+06 cols
damp = 0.00000000000000e+00   calc_var =        0
atol = 1.00e-08                 conlim = 1.00e+08
btol = 1.00e-08               iter_lim =       13
 
   Itn      x[0]       r1norm     r2norm   Compatible    LS      Norm A   Cond A
     0  0.00000e+00   1.286e+02  1.286e+02    1.0e+00  8.0e-03
     1  0.00000e+00   5.589e+01  5.589e+01    4.3e-01  9.4e-01   1.1e+00  1.0e+00
     2  0.00000e+00   3.150e+01  3.150e+01    2.4e-01  3.9e-01   1.8e+00  2.2e+00
     3  0.00000e+00   2.349e+01  2.349e+01    1.8e-01  3.3e-01   2.2e+00  3.4e+00
     4  0.00000e+00   1.624e+01  1.624e+01    1.3e-01  3.2e-01   2.5e+00  5.1e+00
     5  0.00000e+00   1.238e+01  1.238e+01    9.6e-02  2.6e-01   2.9e+00  6.8e+00
     6  0.00000e+00   9.982e+00  9.982e+00    7.8e-02  2.2e-01   3.3e+00  8.6e+00
     7  0.00000e+00   7.004e+00  7.004e+00    5.4e-02  2.4e-01   3.6e+00  1.1e+01
     8  0.00000e+00   5.356e+00  5.356e+00    4.2e-02  1.9e-01   3.9e+00  1.3e+01
     9  0.00000e+00   4.364e+00  4.364e+00    3.4e-02  1.5e-01   4.2e+00  1.4e+01
    10  0.00000e+00   3.322e+00  3.322e+00    2.6e-02  1.8e-01   4.4e+00  1.6e+01
    11  0.00000e+00   2.508e+00  2.508e+00    2.0e-02  1.6e-01   4.6e+00  1.8e+01
    12  0.00000e+00   2.036e+00  2.036e+00    1.6e-02  1.3e-01   4.9e+00  2.0e+01
    13  0.00000e+00   1.607e+00  1.607e+00    1.2e-02  1.5e-01   5.1e+00  2.2e+01
 
LSQR finished
The iteration limit has been reached                      
 
istop =       7   r1norm = 1.6e+00   anorm = 5.1e+00   arnorm = 1.2e+00
itn   =      13   r2norm = 1.6e+00   acond = 2.2e+01   xnorm  = 1.5e+02
In [47]:
fig = plt.figure(figsize=(15,9))
ax1 = plt.subplot2grid((1, 5), (0, 0), colspan=2)
ax2 = plt.subplot2grid((1, 5), (0, 2), colspan=2)
ax3 = plt.subplot2grid((1, 5), (0, 4))

ax1.imshow(g_inv_tot_z1[:, 2].T, cmap='gray', vmin=-5e-1, vmax=5e-1, 
           extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
ax1.set_title(r'$G_{est}$'), ax2.set_xlabel(r'$x_R$'), ax2.set_ylabel(r'$t$')
ax1.axis('tight')
ax1.set_ylim(1.2, 0)
ax2.imshow(g_inv_tot_z1fast[:, 2].T, cmap='gray', vmin=-5e-1, vmax=5e-1, 
           extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
ax2.set_title(r'$G_{est from z0}$'), ax2.set_xlabel(r'$x_R$'), ax2.set_ylabel(r'$t$')
ax2.axis('tight')
ax2.set_ylim(1.2, 0)
ax3.plot(g_inv_tot_z1[nr//2, 2, nt-1:]/g_inv_tot_z1.max(), t, 'r', lw=5)
ax3.plot(g_inv_tot_z1fast[nr//2, 2, nt-1:]/g_inv_tot_z1.max(), t, 'k', lw=3)
ax3.set_ylim(1.2, 0);

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15,9))
ax1.imshow(g_inv_minus_z1[:, 2].T, cmap='gray', vmin=-5e-1, vmax=5e-1, 
           extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
ax1.set_title(r'$G^-_{est}$'), ax2.set_xlabel(r'$x_R$'), ax2.set_ylabel(r'$t$')
ax1.axis('tight')
ax1.set_ylim(1.2, 0)
ax2.imshow(g_inv_minus_z1fast[:, 2].T, cmap='gray', vmin=-5e-1, vmax=5e-1, 
           extent=(r[0,0], r[0,-1], t[-1], -t[-1]))
ax2.set_title(r'$G^-_{est from z0}$'), ax2.set_xlabel(r'$x_R$'), ax2.set_ylabel(r'$t$')
ax2.axis('tight')
ax2.set_ylim(1.2, 0);